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Chaotic dispersal of tidal debris 
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ABSTRACT 

Several long, dynamically cold stellar streams have been observed around 
the Milky Way Galaxy, presumably formed from the tidal disruption of globu¬ 
lar clusters. In integrable potentials—where all orbits are regular—tidal debris 
phase-mixes close to the orbit of the progenitor system. However, the Milky 
Way’s dark matter halo is expected not to be fully integrable; an appreciable 
fraction of orbits will be chaotic. This paper examines the influence of chaos 
on the phase-space morphology of cold tidal streams. Streams even in weakly 
chaotic regions look very different from those in regular regions. We hnd that 
streams can be sensitive to chaos on a much shorter time-scale than any standard 
prediction (from the Lyapunov or frequency-diffusion times). For example, on 
a weakly chaotic orbit with a chaotic timescale predicted to be >1000 orbital 
periods (>1000 Gyr), the resulting stellar stream is, after just a few lO’s of or¬ 
bits, substantially more diffuse than any formed on a nearby but regular orbit. 
We hnd that the enhanced diffusion of the stream stars can be understood by 
looking at the variance in orbital frequencies of orbit ensembles centered around 
the parent (progenitor) orbit. Our results suggest that long, cold streams around 
our Galaxy must exist only on regular (or very nearly regular) orbits; they po¬ 
tentially provide a map of the regular regions of the Milky Way potential. This 
suggests a promising new direction for the use of tidal streams to constrain the 
distribution of dark matter around our Galaxy. 
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1. Introduction 

The dark-matter haloes of galaxies are thought to be triaxial in shape with axis ratios 
and alignments that depend on radius. However, despite suggestive evidenee from a range of 
eomplimentary observational methods, this fundamental predietion from ACDM eosmology 
has not been eonelusively verihed. Around other galaxies, it is generally hard to measure the 
3D mass prohle beeause informative traeers are rare and the haloes are seen in projeetion. 
From the Earth’s position within the Milky Way, our view of our own halo and proximity 
gives us a unique ehanee to direetly measure the 6D positions of stars and model the shape of 
the mass distribution at large radii. The Milky Way halo has a low density of visible traeers, 
but luekily many of the halo stars are likely assoeiated with debris stripped from stellar 
systems and thus eontain additional information eneoded by the formation meehanism. 

As a satellite galaxy or globular eluster orbits within some larger system, mass is eroded 
by tidal forees from the host-galaxy potential. Along regular, mildly eeeentrie orbits, mass 
stripped from the progenitor has a small spread in orbital properties (e.g., energy, angular 
momentum). Onee the debris has evolved far enough from the progenitor system that the 
self-gravity of the progenitor ean be ignored (usually a fast proeess relative to the orbital 
time), the stars evolve essentially as an ensemble of test partiele orbits in the potential of the 
host system. The debris remains eoherent as a tidal stream if the phase-mixing time-seale 
is long: a small ensemble of regular orbits reaehes a fully phase-mixed state in a timeseale 
^ ctq^, where cr^ is the dispersion in fundamental frequeneies of the ensemble (tidal debris 
from a globular eluster typieally has frequeney spreads p^O.1-1%, so it ean take hundreds to 
thousands of orbital periods to fully phase-mix). The ensemble spreads due to shearing from 
slight variations in their fundamental frequeneies, whieh, for tidal streams, preferentially 
oeeurs along one dimension (Merritt & Valluri 1996; Helmi & White 1999). 

The morphologieal (density) evolution of the tidal debris therefore depends on the spread 
of orbital properties (e.g., aetions or frequeneies) of the debris and the orbit of the progenitor 
system, both of whieh are also determined by the shape and radial prohle of the gravitational 
potential of the host galaxy. By modeling the observed phase-spaee density of stream stars 
along with the host galaxy potential it is hoped that we may infer the 3D mass distribution 
of the host. Many tidal streams are observed around the Milky Way, M31, and other nearby 
galaxies; the known streams span a large range of distanees—from ~ 10 to 100 kpe—and 
progenitor masses—from ~10^ to ~10® Mq in stellar mass— (Ibata et al. 1994; Odenkirehen 
et al. 2001; Belokurov et al. 2006; Grillmair & Dionatos 2006; Grillmair 2006; Bonaea et al. 
2012). There has been extensive work on developing methods to use data from these streams 
to measure properties of the Milky Way’s dark matter halo. These methods span a range of 
eomplexity from orbit-htting, to Streakline (Kiipper et al. 2012) or ‘partiele-spray’ models 
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(Gibbons et al. 2014), to action-space density modeling (e.g., Sanders 2014; Bovy 2014), to 
N-body simulations (e.g.. Law & Majewski 2010). All methods have been tested in some 
way on simulated observations of data and these tests typically demonstrate the recovery of 
parameters for analytic, static potential forms. 

One example of stream modeling in a multi-component (static, analytic) potential was 
presented by Pearson et al. (2015), who aimed to reproduce observations of the stellar stream 
density from the globular cluster Palomar 5 in a single oblate and single triaxial potential 
using Streakline (Kiipper et al. 2012) and N-body models. They used the observed number 
density of stream stars, a limited number of radial velocities for stream members, and the 
sky position, distance, radial velocity, and proper motion of the cluster itself to ht model 
streams to the data. In the oblate potential (a three-component bulge+disk+spherical halo 
potential), a thin model stream was easily found that reproduces the observed stellar density 
morphology of the stream. In the triaxial potential (the potential from Law & Majewski 
(2010): a three component bulge+disk+triaxial halo ht to Sagittarius stream data), the 
model streams generically formed large, two-dimensional ‘fans’ of debris near the ends, and 
no physically reasonable progenitor orbits could be found that reproduced the observed 
thinness and curvature of the stream given the observational constraints of the present-day 
position and velocity of the cluster. The result in Pearson et al. (2015) demonstrates that the 
morphology of a stream alone can be used to rule out a potential. With an understanding of 
the circumstances that lead to the differences in stream morphology, this could be a powerful 
tool for rejecting potentials from positional information alone. 

The obvious difference between the two potentials considered by Pearson et al. (2015) 
is the extra symmetry of the oblate potential. It is well known that the number of degrees of 
freedom of a potential plays a critical role in determining the orbit structure of the potential; 
Hamiltonians with more than two degrees of freedom generically contain significant chaotic 
regions. Pearson et al. (2015) tested the stochasticity of the orbit of the progenitor that 
produced ‘fanned’ debris by computing the Lyapunov exponent along this orbit but found 
that it is consistent with being regular over dynamically relevant timescales (many Hubble 
times). It has been shown previously that along some strongly chaotic orbits and in live 
cosmological haloes, tidal streams do form large, diffuse ‘fans’ of debris (e.g., Fardal et al. 
2015; Ngan et al. 2015), however it is unknown how the resultant properties of the debris 
(e.g., density or length of the stream) depend on the degree of stochasticity. The result from 
Pearson et al. (2015) suggests that even weak chaos (as measured by the Lyapunov exponent) 
may affect the density evolution and therefore observability of tidal streams. Understanding 
why this occurs and developing a measure to quantify the importance of this enhanced 
density evolution is a promising new direction to be explored further. 
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In this Article^ we study the effeet of ehaotie diffusion of the fundamental frequeneies of 
individual orbits on the density evolution of tidal debris. We ehoose a simple, eosmologieally 
motivated model for a triaxial potential, analyze the degree of ehaos as eomputed from single¬ 
orbit diagnosties for grids of eonstant-energy orbits, and eompare these results to measures 
of the density evolution of finite-volume ensembles of orbits (meant to mimie tidal debris). 
We find that even when the ehaotie timeseale is predieted to be long for a given orbit, ehaos 
may manifest over mueh shorter times in small orbit ensembles through the ehaotie diffusion 
of the eonst it uent orbits. For a ehaotie orbit, the frequeney speetrum of the orbit evolves 
with time: for a small ensemble, the spread in frequeneies is therefore time-dependent, whieh 
eould enhanee phase-mixing. This idea supports a reevaluation of the importanee of ehaos in 
galaetie haloes and implies that the amount of ehaos in a given potential may have signiffeant 
eonsequenees for the observability and survivability of thin, eold tidal streams. 

This paper is organized as follows. We review relevant nonlinear dynamies in Seetion 2. 
In Seetion 3, we deseribe our ehoiee of potential, method for numerieal orbit integration, and 
introduee the ehaos indieators used in this work. Our results are split into three subseetions: 
in Seetion 4.1 we present iso-energy grids of orbits and diseuss the orbit elasses and ehaotie 
timeseales present in our potential; in Seetion 4.2 we study the density evolution of small 
ensembles of orbits around eaeh orbit of the previous seetion; in Seetion 4.3 we deseribe the 
behavior of ehaotie diffusion and use this to explain how ehaos is relevant for tidal streams 
over short times. We diseuss the implieations of our results in Seetion 5, and eonelude in 
Seetion 6. 


2. Review of nonlinear dynamics 


To explore the question of if and how ehaos manifests in the density evolution of orbit 
ensembles over timeseales mueh shorter than that predieted from generie ehaos indieators, 
we must first understand the behavior of individual orbits in eomplex gravitational potentials 
and the orbital struetures in the potential itself (i.e. the strength of resonanees and ehaos). 
An orbit in an N degree of freedom (dof) Hamiltonian, H, is a set of 2N quasi-periodie time 
series, 

(ti)i(i), W2jv(i)) = (9i(i), qNit),Piit), ■■■,PN{t)) (1) 

where qk and pk are eonjugate eoordinates in the sense that 
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for all t. If bounded, the motion in any eomponent, Wk{t)^ ean be deseribed as a Fourier 
sum, 

oo 

U)k{t) (4) 

3 

where the a^j are eomplex amplitudes. 

A regular orhit is a set of sueh time series that ean be transformed to a speeial set 
of eanonieal eoordinates known as angle-aetion eoordinates. In these eoordinates, the posi¬ 
tion variables are angles, 0, that increase linearly with time with rates set by N constant, 
fundamental frequencies, = (Qi, The frequency of a Fourier component (the ujk 

in Equation 4) for any individual component of motion are just linear, integer combina¬ 
tions of the fundamental frequencies, ft —that is, for a regular orbit, any Fourier component 
frequency may be written 


uij = ■ n (5) 

where is a vector of N integers. The conjugate momentum coordinates—the actions, 
J —are constants of motion. Even stronger, the actions are isolating integrals and any pair 
are in involution such that 

[Ja, Jf<] = 0 (6) 

where [•, •] is the Poisson bracket. This implies that for an N dof system, a regular orbit has N 
independent constants of motion and the motion is therefore restricted to an A^-dimensional 
manifold embedded in the 2N dimensional phase space. The topology of angle-action space 
is toroidal and any regular orbit in an N dof Hamiltonian can be understood as motion on 
the surface of an A^-torus. Each set of actions, (Ji,..., Jn) (or frequencies), labels a torus, 
and regular orbits are sometimes referred to in terms of their orbital tori (see, e.g.. Section 
50 in Arnold 1978, Section 10-5 in Goldstein 1980, Section 3.5 in Binney & Tremaine 2008). 


2.1. Orbits in integrable potentials 

A Hamiltonian or potential is said to be globally integrable when the number of isolating 
integrals of motion is equal to the number of degrees of freedom and a transformation to 
angle-action coordinates may be done globally—for example, the transformation to angle- 
action coordinates may be written as a function of arbitrary phase-space coordinates and 
the functional form is independent of phase-space position (e.g., Goldstein 1980). The con¬ 
dition for global integrability is very restrictive and the likelihood that a Hamiltonian is 
globally integrable decreases as the number of degrees of freedom increase (e.g., Lichtenberg 
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& Lieberman 1983) and as the potential beeomes more eomplex (i.e. eontaining multiple 
eomponents with different shapes). In a globally integrable potential, the Hamiltonian may 
be written solely in terms of the actions, H = H{J). Galactic potentials are almost certainly 
not globally integrable but it is useful to understand the orbit structure in integrable systems 
before extending to more general potentials. For example, there are four general classes of 
orbits in the ‘Perfect Ellipsoid’ potential (an integrable triaxial potential and special case 
of the Stackel potential; see, e.g., Kuzmin 1973; de Zeeuw 1985): box, inner long-axis tube, 
outer long-axis tube, and short-axis tube orbits. Regular orbits in non-integrable triaxial 
potentials can typically still be identiffed with these classes. Tube orbits have a non-zero 
time-averaged angular momentum about either the long or short axis and therefore never 
pass through the center of the potential (hence are centrophobic). Box orbits instead have 
a zero time-averaged angular momentum and therefore have ffnite probability of passing 
through the center of the potential. These orbits are generally centrophilic, though some 
resonant box orbits are also centrophobic (e.g., banana, pretzel, ffsh). 

The frequencies of a generic orbit are typically incommensurable—that is, n Vt ^ 0 
for any integer vector, n, with reasonable magnitude.^ These non-resonant orbits uniformly 
cover the surface of an orbital torus. If instead there exists a relation of the form n-il = 0 the 
orbit is referred to as a resonant orbit. Resonant orbits are confined to a surface with lower 
dimensionality than the surface of an orbital torus, depending on the number of resonance 
relations obeyed; in a triaxial potential, orbits may obey either zero, one, or two resonance 
relations. We refer to orbits that obey a single resonance relation as uni-resonant orbits, 
and if an additional resonance relation exists, bi-resonant. Uni-resonant orbits in a triaxial 
potential are confined to a 2D surface, and bi-resonant orbits are closed ID curves. The 
resonant structure of a potential—the relative importance of particular resonance integer 
vectors—is difficult to compute, but determines the global behavior of orbits in the potential. 
In plots of frequency ratios, the resonances appear as lines; Figure 1 (left panel) shows a 
cartoon portrait of a portion of frequency-space for an integrable potential with example 
resonance lines, non-resonant, and resonant orbits marked. For an integrable potential, all 
orbits are regular. 


more precise condition is stated in terms of a diophantine condition, e.g., |n • i7| > o; |n| where 

a, 7 > 0. 
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Integrable Near-integrable Non-integrable 




Fig. 1.— An illustrative demonstration of the orbit strueture for integrable (left), near- 
integrable (middle), and non-integrable (right) potentials in terms of the ratios of the funda¬ 
mental frequeneies. In an integrable potential, only resonant and non-resonant orbits exist, 
and all orbits are regular (resonanees appear as lines in frequeney-ratio-spaee). If the poten¬ 
tial is perturbed mildly, resonanee layers form around the resonanees that host near-resonant 
orbits whieh behave like resonant orbits but have an additional frequeney eorresponding to 
libration about the resonanee. Where resonanees overlap or near separatriees, stoehastie 
layers form and orbits will be ehaotie. For more strongly perturbed potentials, many of the 
non-resonant orbits may beeome ehaotie but resonanee layers may grow and still remain 
regular. 


2.2. Orbits in near- and non-integrable potentials 

The orbit strueture of near-integrable potentials ean be understood by eonsidering a 
Hamiltonian that is a small perturbation away from being globally integrable—that is, a 
Hamiltonian that may be written 

H(J,d) = Ho{J)+eH,(J,0) (7) 

where Hq{J) represents an integrable Hamiltonian, and e is a small parameter that deter¬ 
mines the perturbation strength (a deseription of perturbation theory applied to nonlinear 
Hamiltonians is given in Liehtenberg & Lieberman 1983). When 0 < |e| ^ 1, resonant sur- 
faees beeome ‘thiek’ resonant layers, within whieh orbits are qualitatively similar to the par¬ 
ent resonant orbit (e.g., Merritt & Valluri 1999). These resonant layers are then generieally 
surrounded by stoehastie layers where ehaotie motion oeeurs (in the vieinity of separatriees). 
Chaotie orbits are not bound to the surfaee of a torus and instead hll a (2N-l)-dimensional 






volume where a eontinuous set of tori would exist in an unperturbed Hamiltonian, Hq. While 
the frequeney speetrum of sub-seetions of a regular orbit are indistinguishable (i.e. measured 
over a hnite interval of time), the frequency spectra of sub-sections of a chaotic orbit will 
evolve stochastically. 

Figure 1 (middle panel) shows a cartoon of frequency space for a near-integrable poten¬ 
tial (a small perturbation away from the potential in the left panel) assuming the resonances 
are stable. Much of the structure that was present in the integrable potential remains in 
the near-integrable case, but the differences are highlighted. Orbits in the resonant layers 
surrounding the resonances (grey) are near-resonant orbits that librate around the resonance 
and have hnite thickness (e.g., Merritt & Valluri 1999). Chaotic orbits in the stochastic layer 
(red) behave erratically depending on the surrounding resonance structure. If the stochastic 
layer is small and the chaotic orbit is therefore conhned, the orbit may behave nearly regular 
for long periods of time. If the resonance in the unperturbed Hamiltonian is unstable, all 
orbits associated with the resonance will be chaotic; unstable resonances form linear gaps in 
frequency-space (this is not shown in the cartoon but are seen in Figure B.l). Thus, only 
some resonances in the perturbed system will retain the signature of a resonance. 

For small values of e, many regular orbits survive and only small chaotic regions are 
introduced, especially in the vicinity of the intersection of two resonance lines (commonly 
referred to as ‘resonance overlap’; see Chirikov 1960). As the strength of the perturbation in¬ 
creases, eventually most tori associated with non-resonant motion will be destroyed. Figure 1 
(right panel) qualitatively shows this phenomenon—the resonant and resonant-layer orbits 
may still be regular, but many or all of the non-resonant orbits will be chaotic. As the pertur¬ 
bation strength increases, eventually the uni- and bi-resonant tori are also destroyed—these 
are less susceptible to destruction from perturbations (for a more quantitative illustration of 
this transition from integrability to global chaos, see Figure 9 in Valluri & Merritt 1998). 

When e is large, there is no general prediction for the resulting behavior, however it 
seems that more complicated and physically motivated triaxial potential models for galaxies 
follow the intuition gained from the small-perturbation picture described above, at least for 
certain parameter choices (e.g., Valluri & Merritt 1998; Merritt & Valluri 1999). We therefore 
expect a large number of regular orbits will survive—the so-called KolmogorovArnoldMoser 
(KAM) tori—however the tori that survive will be separated by regions of chaotic motion. 
Any transformations to angle-action coordinates must be dehned local to each resonance 
region due to the destruction of tori and chaotic motion which lead to discontinuous changes 
in orbital properties. 
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2.3. The behavior of chaotic orbits in non-integrable potentials 

Chaotic orbits have no orbital actions and only conserve energy (if the potential is 
time-independent). The orbits therefore do not have a single set of fundamental frequencies, 
but rather the frequencies that describe the character of motion evolve with time. In near- 
integrable potentials, the frequencies of consecutive sub-sections of a chaotic orbit diffuse 
both around resonance layers^ For weakly chaotic orbits, motion around a resonance layer 
can occur with a frequency close to the libration frequency of the nearby stable orbits in 
the resonance layer. Thus, if the resonance libration frequency is small, motion around a 
resonance can modulate the frequency spectrum of an orbit over an orbital time. However, 
the stochastic layers are often bounded in the direction orthogonal to the resonance by other 
stable, resonant regions so that the frequencies or actions can not change by large factors 
(unless there are other nearby resonances and overlapping stochastic layers, in which case 
the motion may be strongly chaotic). 

The rate of diffusion along stochastic layers via Arnold diffusion depends on the local 
resonant structure and is hard to predict. This has been done analytically for simple poten¬ 
tials (e.g., Chirikov 1979). For systems with N > 2 dof, the stochastic layers connect and 
form an intricate network of stochasticity known as the Arnold web; an orbit that ergodically 
mixes over its energy hypersurface must traverse this web, though the timescales typical for 
this phenomenon are many thousands of orbital periods. 

Arnold diffusion is not expected to be significant for most orbits over timescales relevant 
to galaxies (10s of orbits), however chaotic motion across resonances can occur over short 
times. If a stochastic trajectory is surrounded by regular orbits, it may be trapped around a 
regular parent orbit for long periods of time before escaping to another such semi-bounded 
region where it can become trapped around another parent orbit (a process by which it may 
eventually explore the whole Arnold web)—such orbits are commonly referred to as ‘sticky 
orbits.’ Additionally, if the volume of the surrounded region in frequency space is comparable 
to the characteristic spread of frequencies in the tidal debris, this small-scale evolution will 
be important for tidal debris. 


^Note that during this diffusion, the frequencies never exactly hit those of a KAM torus but evolve 
stochastically around these discrete, stable tori (cf. Figure 2 in Laskar 1999). (a sort of stochastic libration) 
and along the stochastic layers that surround the resonances (Arnold diffusion). 
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2.4. Mixing of orbit ensembles 

An ensemble of regular orbits (e.g., tidal debris) will phase-mix beeause of (small) 
differenees in the fundamental frequeneies of the orbits. The frequeney distributions of thin 
streams are generally elose to one-dimensional—that is, one eigenvalue of the distribution 
of frequeneies for an orbit ensemble will be mueh larger than the others beeause of the 
loeal shape of the Hessian of the potential (Helmi & White 1999; Sanders & Binney 2013; 
Bovy 2014). Phase-mixing generieally leads to power-law deeay of the mean density of the 
ensembles: initially, the density deereases linearly in time beeause of the large, nearly one¬ 
dimensional spread in frequeneies, then may proeeed as to depending on relative 
sizes of the other eigenvalues of the frequeney-spaee distribution (e.g., Helmi & White 1999; 
Vogelsberger et al. 2008). 

Generieally, a small ensemble of ehaotie orbits will lose eoherenee mueh faster than for 
regular orbits (see, e.g., Kandrup & Mahon 1994; Merritt & Valluri 1996; Kandrup & Siopis 
2003) , however this depends on the details of the resonant strueture around the ensemble and 
the ehaotie evolution of the individual orbits and is thus diffieult to prediet. For example, 
ensembles of orbits ‘stuek’ between resonanees may quiekly spread to hll the allowed volume, 
but then the orbits must eseape this eonhnement and diffuse through the Arnold web to reaeh 
a fully mixed state (Merritt & Valluri 1996). That is, while the orbit is stuek, the small-seale 
variations effeetively eause an inerease in the varianee of the frequeney distribution, whieh 
would enhanee mixing of the debris in eonffguration-spaee. In this work, we investigate the 
eonsequenees of short-time but small-seale frequeney evolution and hypothesize that this may 
explain the enhaneed density evolution of tidal debris around weakly ehaotie regions where 
ehaotie timeseales are predieted to be long (e.g., Pearson et al. 2015). We then diseuss how 
this would affeet our understanding of the eoherenee and density evolution of tidal streams. 


3. Numerical methods 

Our goal is to map the orbit strueture of arbitrary (galaetie) potentials, with an emphasis 
on identifying the ehaotie orbits and understanding the evolution of these ensembles of orbits 
over short times. In partieular, we aim to understand how this ehaos-enhaneed density 
evolution ean affeet tidal stream morphology. In this seetion, we deseribe the methods we 
will use to deteet and quantify the strength of ehaos for large grids of orbits. 
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3.1. Potential choice 

The density distributions within dark-matter haloes formed in eosmologieal N-body 
simulations are generieally triaxial (e.g., Jing & Suto 2002; Bett et al. 2007; Zemp et al. 2009; 
Vera-Ciro et al. 2011). With the inelusion of baryonie physies and sub-grid preseriptions 
for energy input due to supernovae and other feedbaek meehanisms, the inner potential 
(<0.1i?vir for a p^lO^^ Mq halo mass) typieally beeomes more spherieal, though the magnitude 
of this reshaping depends on the partieular merger history and star formation effieieney 
within a given halo and the mass and shape of the baryonie eomponent (e.g., Dubinski 1994; 
Kazantzidis et al. 2004; Debattista et al. 2008; Bryan et al. 2013; Butsky et al. 2015, though 
in Milky Way-like galaxies, baryonie disks will add non-spherieity to the total potential). It 
is less elear what happens to the outer halo (e.g., Zemp et al. 2011; Valluri et al. 2013). 

Jing & Suto (2002, hereafter JS02) found that a triaxial generalization of the NFW 
density prohle (Navarro et al. 1996) generates exeellent hts to the density distributions within 
haloes in their high-resolution (dark-matter-only) N-body simulations, and they provide 
probability distributions for the axis ratios of a large sample of these haloes. JS02 hnd median 
axis ratios of c/a ^ 0.55 and b/a ^ 0.77 where a is the major axis, b the intermediate, and c 
the minor axis.^ These are largely eonsistent with hndings from more reeent simulations (e.g., 
Bett et al. 2007; Vera-Ciro et al. 2011; Butsky et al. 2015; Zhu et al. 2015) and eonsistent 
with eonstraints from weak leasing that plaee a lower limit on minor-to-major axis ratios 
of c/a > 0.5 (van Uitert et al. 2012). JS02 hnd signiheant seatter in the distributions of 
eoneentration parameter, Cg, or seale radius (depending on ehoiee of parametrization). 

All of these parameters are speeihed in terms of the density; for orbit analysis, we 
need to determine the form of the potential in terms of these parameters, whieh, in general, 
requires numerieal integration of the density at eaeh position of interest. For eomputational 
effieieney, many authors instead express the triaxiality in the form of the potential, but this 
ean lead to unphysieal situations where the density beeomes negative. Lee & Suto (2003) 
derive a perturbative expansion of the potential integral for a triaxial NFW density and show 
that the expansion is aeeurate even for modest axis ratios (e.g., the median values shown 
above). 

In this work, we use the triaxial potential expression from Lee & Suto (2003), parametrized 
in a slightly different manner. In terms of spherieal eoordinates^ with the radius normalized 


^Note that JS02 use the opposite notation so that c is the major and a is the minor axis, 
^(r, b, 6*) = (radius, azimuth, colatitude) 
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Fig. 2.— Equipotential contours (left) and isodensity contours (right) for the triaxial NFW 
potential considered in this work. For the potential plot there are eight contour levels evenly 
spaced and linear in the value of the potential. For the density plot there are eight contour 
levels logarithmically spaced from 10^ Mq kpc“^ to 10^ Mq kpc“^. 


by the scale radius, u = r/vs 


$(n, (/), ^) 

A = 


^ Fi(«) + ^{el + el)F2{u) + ^[(e6sm6'sm<j!>)2 + {ecCosd)^]Fi{u) 

(in2 - 2) + (l“^ ~ 4) ^ 


( 8 ) 

(9) 


where — {b/ay, Cc = y/l — {c/a)^, and Vc is the circular velocity at the scale radius, 

r^, for the spherical case. The functions Fi{u) are given in the appendix of Lee & Suto (2003). 
We chose = 20 kpc and Vc = 175 km s“^ by taking the mean halo concentration for a 
Myir ~ 10^^ M 0 halo, Ce ~ 5, from Jing & Suto (2002) and by assuming Ryir ps 200 kpc. 
Figure 2 shows equipotential contours of this potential, and Table 1 summarizes the potential 
parameters. 


This potential is a simple (and unrealistic) model for the total potential of a Milky-Way- 
like galaxy, however it represents a conservative choice for exploring the structure of orbits 
in the haloes of such galaxies. Realistic galactic potentials will have a signihcant component 
due to the disk and bulge, radially changing axis ratios or orientations (e.g., Romanowsky 
& Kochanek 1998; Kazantzidis et al. 2004; Debattista et al. 2008; Vera-Ciro et al. 2011; 
Butsky et al. 2015), signihcant substructure (Moore et al. 1998; Zemp et al. 2009), or time 
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Parameter 

Value 

Vc 

175 km s-i 

Ts 

20 kpc 

a 

1 

b 

0.77 

c 

0.55 

Tabc 

0.58 


Table 1: Summary of parameters for the triaxial NFW potential (Equation 8) used in this 
work. Triaxiality is introduced in the density (rather than the potential) to ensure that the 
density is physical at all radii. Velocity scale, scale radius, and axis ratios are chosen to 
match the median halo parameters for a Mvir ~ 10^^ Mq halo from (Jing & Suto 2002). The 
triaxiality parameter, Tabc = is also given. 


dependence (either from bulk rotation, mass growth, mergers, etc.; see, e.g., Bailin et al. 
2005) . We expect inclusion of any of these effects to increase the complexity of the resonant 
structure and inffuence of chaos (see Section 5). 


3.2. Orbit integration 

We use the Dormand-Prince 8th-order Runge-Kutta scheme (Prince & Dormand 1981) 
to integrate orbits in the above potential. Speciffcally, we use a Python wrapper over the C 
implementation by Hairer et al. (1993). For all orbits we ensure that energy is conserved to 
\AE/Eo\ < 10“® by the end of integration, however most orbits conserve energy to a part 
in ^10“^^. Unless otherwise speciffed the integration timesteps are chosen so that there are 
512 steps per strongest orbital period component, but the integrator uses adaptive stepping 
between each main step in order to satisfy a speciffed tolerance (we set the absolute tolerance 
to ^100 times machine precision, atol = 10“^^). 


3.3. Lyapunov exponents 

The most well-known method for assessing chaotic motion is to analyze the Lyapunov 
spectrum or maximum Lyapunov exponent (MLE) of an orbit (Lyapunov 1992). The MLE 
measures the mean rate of divergence of two inffnitesimally separated orbits and is only 






-14- 


strictly defined in terms of a limit that goes to infinite time. Thus, we ean never truly 
eompute the MLE and it ean take integration for many thousands of orbital periods to 
compute a converged numerical approximation of the MLE for a moderately chaotic orbit. 
In this work, we use the algorithm introduced by Wolf et al. (1985) for computing the MLE 
(for more a more detailed description of this algorithm, see Appendix A). 

The MLE, A^ax, is interpreted as a rate that quantifies the exponential divergence of 
infinitesimally close chaotic orbits. It is therefore useful to consider the corresponding e- 
folding time by inverting the rate, 

( 10 ) 

^max 

We will use this as the prediction from the Lyapunov exponent for the timescale over which 
chaos should be dynamically important for a given orbit. 


3.4. Frequency diffusion rate 

Bounded, regular orbits in a triaxial potential have three fundamental frequencies, fi, 
that determine the periodic behavior of motion. The motion in any canonical coordinate 
can therefore be decomposed as a Eourier sum (Equation 4) where the Eourier frequencies 
are linear, integer combinations of the fundamental frequencies (Equation 5). Laskar (1993) 
introduced a method for recovering the fundamental frequencies of an orbit that effectively 
uses fast-Eourier transforms (EETs) of complex combinations of the motion (e.g., x{t) + 
iVx{t)) to identify the frequencies. This method is referred to as ‘Numerical Approximation 
of Eundamental Erequencies’ (NAEE) and has been used extensively in planetary dynamics 
(e.g., Laskar & Robutel 1993; Laskar 1996) and galaxy dynamics (Papaphilippou & Laskar 
1998; Valluri & Merritt 1998), especially in the study of orbits in triaxial systems. We have 
implemented and tested a version of this procedure in the Python programming language. 
Our implementation differs slightly from the original definition and from that used in Valluri 
& Merritt (1998); we refer to this slight modification of the algorithm as SuperFreq (Price- 
Whelan 2015) and the code is open-source and publicly available on GitHub.^ Eor more 
details about the algorithm and differences with previous work, see Laskar (1988, 1993); 
Papaphilippou & Laskar (1996) and Appendix B. 

If an orbit is chaotic, the motion can no longer be expanded in terms of a single set 
of fundamental frequencies. Eor a weakly chaotic orbit, the orbit may appear consistently 
periodic over long windows of time. SuperFreq will pick out a set of frequencies for chaotic 


^https://github.com/adrn/SuperFreq 
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orbits that correspond to the largest peaks in the power spectrum of the orbits, however 
these peaks will change location and amplitude with time. For more strongly chaotic orbits, 
the power spectrum will be quite noisy and the peak frequencies may change erratically 
when comparing two consecutive sections of orbit. The frequencies picked out by SuperFreq 
for such orbits will therefore represent the average periodic nature of the orbit over a given 
integration window. Following previous work, we dehne the fractional frequency diffusion 
rate, IZ, in the kth fundamental frequency as 




( 1 ) 


At 


( 11 ) 


where the upper index refers to the two consecutive sections of orbit and At is the length 
of each integration window (Laskar 1993; Valluri & Merritt 1998; Valluri et al. 2012). By 
inverting this rate, we can compute the timescale over which we expect order-unity changes 
to the fundamental frequencies: the frequency diffusion time is defined as 


tn = (max IZk) ^ (12) 

ak 

where the maximum is taken with respect to the corresponding amplitudes, a^, of the fun¬ 
damental frequency components (see Valluri et al. 2012). 


For a small ensemble of orbits (e.g., tidal debris), a more relevant timescale is the 
time over which the change in frequencies for a single orbit is comparable to the spread of 
frequencies in the ensemble. We can estimate this timescale by multiplying the frequency 
diffusion time by a factor equal to the fractional spread in frequencies of the debris. For 
example, a globular cluster typically has ^0.1% spreads in fundamental frequencies, so by 
multiplying the frequency diffusion time by / = 10“^, we can estimate the time (in number 
of orbital periods) over which we expect the frequencies to evolve by this amount. 


4. Results 

In Section 4.1, we generate grids of orbits in the potential described above to map the 
orbital structure of the potential. We classify each orbit in terms of the strength of chaos 
along the orbit as computed using the Lyapnov and frequency diffusion times of Section 3. 
With this initial classiffcation, in Section 4.2 we follow the evolution of ensembles of trajec¬ 
tories generated around each orbit in the initial grid. We ffnd that the conffguration-space 
density of ensembles around weakly chaotic orbits evolve faster (e.g., in mean density) than 
expected given the timescales over which chaos is computed to be relevant for the parent 
orbits. In Section 4.3, we explain this phenomenon in the context of how chaotic diffusion 
occurs (e.g.. Section 2.3). 
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4.1. Part I: Lyapunov and frequency diffusion times 

We generate an isoenergy grid of initial eonditions along the xz {y = 0) plane® with 
energy (per unit mass), E, ehosen to span a range of distanees eomparable to the seale radius 
of the potential (E = —(397.2 km s“^)^ in physieal units; see Table 1). We hx Vx = Vz = 0, 
and eompute Vy from the energy. This grid generates all of the major orbit elasses (short- 
axis tubes, inner long-axis tubes, outer long-axis tubes, stoehastie intermediate-axis, and 
box orbits). The most numerous orbits in this grid are the short-axis and long-axis tubes 
that eireulate about either the minor or major axis. Thin tidal streams may preferentially 
form along tube orbits rather than box orbits beeause of the faster disruption and debris 
diffusion expeeted for stellar systems on radially plunging orbits. Figure 3 shows the grid of 
initial eonditions in the xz plane—eaeh pixel in this grid represents an orbit, and the pixels 
are eolored by the median perieentrie distanee (left) and median apoeentrie distanee (right) 
over an integration time of 64 orbital periods. From here onwards, all lengths are given in 
units of seale radii, r^, veloeities in units of eireular veloeity, Uc, (see Table 1) and times in 
units of orbital periods, T. 

We integrate all orbits in the grid for 10000 orbital periods and use the method deseribed 
in Seetion 3.3 to eompute the MLEs. Figure 4 again shows the grid of initial eonditions, 
but now the eolor eorresponds to the logarithm of the inverse of the MLE (the Lyapunov 
time). The darker pixels have shorter Lyapunov times and are more ehaotie. Beeause of the 
ffxed the integration time, the MLE eannot deteet weak ehaos and the majority of orbits 
appear to be regular beeause they have exeeedingly long Lyapunov times (all white points 
have (<a/T) > 1000). 

Eor eaeh orbit, we also separately integrate for 256 orbital periods and use SuperFreq 
to eompute the fundamental frequeneies for the two eonseeutive seetions of 128 orbital pe¬ 
riods. We have ehosen this window size so that we reeover the frequeneies for regular orbits 
with fraetional error ~ 10“® (we estimate the error in frequeney reeovery using the method 
deseribed in Laskar 1993). With this integration window we are able to sueeessfully reeover 
frequeneies for >99.9% of the orbits SuperFreq. Eigure 5 shows the same grid of initial 
eonditions as in Eigure 4, but now the greyseale intensity is set by the logarithm of the 
frequeney diffusion time. The darker pixels have shorter frequeney diffusion times and are 
more ehaotie. This map reveals the interseetion of this partieular energy hypersurfaee with 
the rieh strueture of resonant surfaees present in this potential and highlights the aeeuraey 
of SuperFreq-weak ehaos (grey) is deteetable over mueh shorter integration periods using 
frequeney analysis, eompared to the many tens of thousands of orbits it would take to deteet 


x is the major and z the minor axis. 
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Pericenter, r^/r 
0.5 


Apocenter, r^/r^ 
1.5 2.0 



1.0 1.5 2.0 


Fig. 3.— A grid of isoenergy orbits initialized on the xz plane. All distances are normalized 
by the potential scale radius. Each pixel in these panels represents a single orbit, and 
the shading of each pixel corresponds to the median pericentric distance (left) or median 
apocentric distance (right) computed over 64 orbital periods. The central black band in the 
left panel and white band in the right panel are tube orbits with apocenter-to-pericenter 
ratios close to one. The four arrows are explained in Section 4.2 and the caption of Figure 4. 


such features with the maximum Lyapunov exponent. The tube orbits in this potential are 
mostly regular or only mildly chaotic—the largest regular regions are associated with the 
short-axis and long-axis tube orbits—however islands of stronger chaos do appear, especially 
at the intersections of resonances where resonance overlap occurs. 

The strongest chaotic regions (black) appear in both of the above grids (Figure 4 and 
Figure 5). Some of the weakly chaotic unstable resonances do appear in the Lyapunov time 
map—for example, near (xo,^o)/'^s ~ (1.5, 0.6) and (xo,^o)/^s ~ (0.6, 0.4) where there is a 
slight hint of weak chaos (grey) in Figure 4. The details of the resonant structure is revealed 
in the frequency diffusion map from integrations of just 256 orbital periods. While there is 
rich structure and a signiffcant number of weakly chaotic orbits, the majority of the orbits 
have estimated chaotic timescales corresponding to thousands of orbital periods and are 
thus not expected to be relevant for tidal stream evolution. In the next section, we analyze 
the density evolution of ffnite-volume ensembles of orbits around each orbit in the above 
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Fig. 4.— The same grid of orbits as shown in Figure 3, but now eaeh pixel is eolored by 
the logarithm of the Lyapunov time (in units of orbital periods). Orbits are integrated for 
a total of 10000 orbital periods. Chaotie orbits with tx/T > 700 appear regular beeause the 
integration time for eaeh orbit is insuffieient to resolve weak ehaos. Four orbits are pointed 
out with arrows—from top to bottom, these are the strongly ehaotie (D), near-resonant (A), 
weakly ehaotie (C), and non-resonant (B) orbits of Table 2 (see Seetion 4.2). 
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T 



Fig. 5.— The same grid of orbits as shown in Figure 3, but now eaeh pixel is eolored by the 
logarithm of the frequeney diffusion time (again in units of orbital periods). The frequeneies 
are eomputed in two eonseeutive windows, eaeh of whieh has length equal to 40 orbital 
periods. The frequeneies are measured preeisely so that small ehanges in the frequeneies ean 
be deteeted over just p^IOs of orbits. Four orbits are pointed out with arrows—from top 
to bottom, these are the strongly ehaotie (D), near-resonant (A), weakly ehaotie (C), and 
non-resonant (B) orbits of Table 2 (see Seetion 4.2). 
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grids in order to compare the effect of ordinary phase-mixing of tidal debris with potentially 
enhanced mixing due to chaos. We then compare the density evolution of the ensembles to 
the single-orbit chaos indicators computed in this section. 


4.2. Part II: Ensemble properties and mixing 

The Lyapunov and frequency diffusion times measure the timescales over which chaos 
is relevant for a given orbit—that is, these quantities are measures of how inffnitesimal 
deviations will diverge on average from some parent orbit, or of how long it takes for the 
frequencies of a single orbit to change by some amount. Tidal debris is disrupted from 
progenitor systems with ffnite spreads in orbital properties (e.g., energy). For a disrupting, 
globular-cluster-scale progenitor, the typical energy or frequency-space dispersion of the 
debris is 0.1-1% of the progenitor orbital energy (assuming masses of 10^-10® M©; Johnston 
1998), but for a dwarf-galaxy-scale progenitor, the dispersion can be ^10%. In this section, 
we ask whether the Lyapunov or frequency diffusion time predict the timescale over which 
a ffnite phase-space volume (e.g., tidal debris) stays coherent. 

It is computationally intractable to run full N-body simulations for the large grid of 
orbital initial conditions of the previous section and we therefore take a simpliffed approach 
for studying how ffnite-volume debris spreads along each of these orbits. We instead consider 
small ensembles of particles meant to represent debris disrupted from a single tidal disruption 
event. For a given set of orbital initial conditions—the ‘parent’ orbit—we integrate for 128 
orbital periods, ffnd the phase-space position of the minimum pericenter over this time, 
initialize a small ensemble of test particle orbits around this position, and integrate the 
orbits of all test particles for some integration time. We are interested in the degree to 
which chaos enhances the mixing rate of orbit ensembles and we therefore want to isolate 
out the effect of ordinary phase-mixing along regular orbits. We set the physical scale of the 
ensemble by the tidal radius in position and the velocity scale in velocity; the initial spread 
in fundamental frequencies will scale as (m/M)^/^ like the tidal radius and velocity scale 
(e.g., Johnston 1998; Price-Whelan et al. 2014). If (ccq^'^^o) are a set of initial conditions at 
pericenter for a parent orbit, then 5xk,i and 5vk^i are the kth components (/c = 1, 2, 3) of the 
position and velocity deviation vectors for the ffh particle and the magnitude of the offsets 
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are assumed to be Normally distributed away from the parent orbit: 


^^k,i ^ A/” (0, Ttide/V^) 

(13) 

^Vk,i ^ A/'(0, ay/Vs) 

(14) 

/ m \ 

(15) 

IIccoId) 

(m« iwd) 

(16) 


where M(< r) is the mass enelosed of the host potential within radius r, m is the mass 
seale of the ‘progenitor,’ and || • || is the Euelidean norm. We take m = 10^ M© to represent 
globular-eluster-like progenitors, and use the spherieally-averaged enelosed mass of the host 
potential to estimate the above debris seales. 

We start by eonsidering four partieular orbits ehosen from the orbit grid of Seetion 4.1: 
a regular, near-resonant orbit (A), a regular, non-resonant orbit (B), a weakly ehaotie orbit 
(C), and a more strongly ehaotie orbit (D). The orbits were ehosen to be elose on the orbit 
grid so that their orbital properties (e.g., apoeenter, perieenter) are similar, but have different 
frequeney diffusion times; all four orbits are long-axis tube orbits. Figure 6 shows the orbits 
in projeetion over an integration period of 1024 orbital periods. The initial eonditions and 
ehaos diagnosties for eaeh orbit are listed in Table 2. Figure 7 shows hnal positions of test- 
partiele ensembles initialized around the four orbits deseribed above and integrated for 64 
orbital periods. A thin stream forms on the near-resonant orbit (left eolumn), a thin—but 
more two-dimensional—stream forms on the non-resonant orbit (middle-left), a more diffuse 
stream forms on the weakly ehaotie orbit with a slightly ‘fanned’ morphology (middle-right), 
and a two-dimensional, ‘fanned’ stream on the more strongly ehaotie orbit (right). Given 
the long Lyapunov and frequeney diffusion times of the weakly ehaotie orbit (tx/T > 900, 
tci/T ^ 2 X 10^), it is surprising that the density evolution of the ensemble on this orbit 
appears to be more diffuse than the regular orbit stream. 

We verify that these orbit ensembles eapture the nature of more realistie stream for¬ 
mation by running N-body simulations of globular-eluster-mass progenitor systems on these 
same four orbits. We use the Self-Consistent Field (SCF) basis funetion expansion eode 
(Hernquist & Ostriker 1992) to run the simulations, whieh we set up to run from apoeenter- 
to-apoeenter (rather than perieenter-to-apoeenter as in the ensembles), but finish with the 
progenitor in the same loeation as the parent orbit in the ensemble evolution deseribed above. 
Figure 8 shows the ffnal partiele distributions rotated so that the angular momentum of the 
progenitor orbit is aligned with the z-axis. From a eomparison of Figures 7 and 8, it is elear 
that the morphology of the ensembles is visually similar to the ‘oldest’ (first stripped) debris 
in the N-body simulations. 
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Fig. 6.— Four representative orbits ehosen from the orbit grid (e.g., Figure 3): a regular, 
near-resonant orbit (A), a regular, non-resonant orbit (B), a weakly ehaotic orbit (C), and a 
more strongly ehaotic orbit (D). All orbits are long-axis tubes with similar pericenters and 
apocenters. Orbits in this Figure were integrated for 1024 orbital periods and are shown in 
projection. 


Visual inspection of Figures 7 and 8 suggests that chaotic mixing of small orbit ensembles 
affects the conffguration-space evolution of an ensemble over short times, even when the 
predicted chaotic timescales (from the Lyapunov and frequency diffusion times) are long: 
The mean, single-orbit chaos indicators are not well-suited for determining the importance 
of chaotic diffusion on tidal stream evolution. As a quantitative measure of this enhanced 
density evolution, we compute the evolution of the mean conffguration-space density for 
orbit ensembles evolved around the four parent orbits (A,B,C,D) described above. At each 
time step, we use kernel density estimation (KDE)^ with the ensemble of particle positions 
to estimate the conffguration-space density ffeld. We use an Epanechnikov kernel with an 
adaptive bandwidth: at each evaluation of the density, we use 10-fold cross-validation to ffnd 
the optimal kernel bandwidth. We evaluate the density at the positions of each particle, p^, 
and compare to the mean initial density, (po). 


^We use an implementation from the Python package scikit-learn (Pedregosa et ah 2011). 
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Orbit ensembles 



Fig. 7.— Final particle positions after integrating unbound, globular-cluster-sized ensembles 
of orbits generated around each of the three N-body orbits (Figure 8). The ensembles each 
contain 1024 particles and are initialized at pericenter. The particle positions qualitatively 
match the morphology of the ‘oldest’ debris from the corresponding N-body simulations 
(Figure 8). 


Figure 9 shows the mean density of the ensemble particles computed at 256 evenly- 
spaced intervals from the initial ensemble distribution to the distribution after 64 orbital 
periods. Over-laid on each panel are qualitative power-laws (i.e. not ht to the density 
evolution) that demonstrate the expected trends: After initial decay, the mean density 
of the near-resonant orbit (A) should follow whereas the non-resonant orbit (B) will 
transition from t~‘^ to after long times (and thus currently display an intermediate power- 
law index). Given the extremely long chaotic timescale for the weakly chaotic orbit (C), we 
would expect the mean density to follow a simple power-law over long times, however at 
?:il6 orbital periods the density clearly diverges and begins to follow a decaying exponential. 
The mean density along the strongly chaotic orbit (D) evolves stochastically but generally 
follows a decaying exponential. 

Motivated by the noticeable discrepancy between the hnal mean density between the 

















-24- 


ID 

xjVs 

y/^s 

zjvs 

VxIVc 

VyjVc 


n/T 

ta/T 

A 

0.85 

0 

1.303 

0 

0.721 

0 

> 700 

> 10^ 

B 

0.85 

0 

1.152 

0 

0.849 

0 

> 700 

> 10^ 

C 

0.85 

0 

1.27 

0 

0.752 

0 

> 700 

fs 3 X 10® 

D 

0.85 

0 

1.434 

0 

0.597 

0 

8.14 

2.5 X 10'* 


Table 2: Orbital initial conditions from the xz grid for orbits with a range of chaotic 
timescales—A is a regular, near-resonant orbit, B a regular, non-resonant orbit, C a weakly 
chaotic orbit, and D a strongly chaotic orbit. Positions (x, y, z) are given in units of scale 
radii, r^, and velocities (v^, Vy, v^) in units of scale velocity, Vc- Chaotic timescales are 
expressed {tx, to) in number of orbital periods. Recall that the frequency diffusion time, to, 
is the time over which we expect order-unity changes in the fundamental frequencies, hence 
why the timescales appear quite long. 


two regular (A,B) and the weakly chaotic (C) orbit ensembles—even though the chaotic 
timescale of the weakly chaotic ensemble parent orbit is several times the integration period 
used above—we compute the hnal mean density for orbit ensembles generated around each 
orbit in the grid described in Section 4.1. For all ensembles, we integrate the orbits for 64 
(parent) orbital periods and compute the initial and hnal values of the mean, conhguration- 
space density. Figure 10 again shows the grid of initial conditions (Section 4.1, but now the 
greyscale indicates the ratio of the hnal mean density to the initial density. Interestingly, 
much of the structure that is visible in the upper-half of the frequency dihusion time map 
(Figure 5) is again visible in this map of the density evolution of orbit ensembles. Many 
features appear as lighter, curved features in the upper portion of this grid where the density 
remains systematically higher—these are the stable resonances of the potential. Darker 
regions have systematically lower densities and correspond to regions of resonance overlap 
where weakly chaotic orbits are found. 
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iV-body simulations 



Fig. 8.— We use the Self-Consistent Field (SCF) basis function expansion code (Hernquist 
& Ostriker 1992) to run N-body simulations of globular-cluster-mass progenitor systems 
on the three orbits of Figure 7. The progenitor in each simulation is initialized as a 10^ 
particle Plummer sphere and the background triaxial NFW potential is turned on slowly 
over 250 Myr to reduce artihcial gravitational shocking. We start the progenitor systems 
at apocenter and integrate for ^64 orbital periods so that each simulation hnishes again at 
subsequent apocenter. The mass of the progenitor is set to 10^ Mq, and the length-scale 
of the Plummer sphere is set to 5 pc to get ^50-75% mass loss over the integration time. 
Panels show particle positions from the hnal snapshots of the simulations—for visualization 
the positions have been rotated so that the angular momentum of the progenitor at the hnal 
snapshot are aligned with the z-axis. These simulations conhrm that the test-particle orbit 
ensembles (Figure 7) do (qualitatively) capture the nature of the early-stripped tidal debris. 
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Fig. 9.— Evolution of the mean density (normalized by the initial density) of orbit ensembles 
around the four orbits (Figure 6) over 64 orbital periods. Over-plotted are various power- 
laws that qualitatively show the deeay of the mean density for the ensembles: the mean 
density of the near-resonant orbit (A) follows nearly the non-resonant orbit (B) is 
slowly transitioning from to the weakly ehaotie orbit (C) initially follows but 
then exponential deeay takes over; the strongly ehaotie orbit (D) is roughly exponential with 
fairly erratie density evolution. 


However, it is surprising that these features are present: The ehaotie timeseales predieted 
from both the Lyapunov and frequeney diffusion times were typieally 100 to 1000s of orbital 
periods for many of the unstable features around resonanees. The timeseales in the lower 
portion of the grid (where the orbits are primarily stable, short-axis tube orbits) are simply 
too long to deteet density differenees over 64 orbital periods even from weak ehaos in this 
partieular ehoice of potential. 

It is worth noting that the frequeney diffusion time map prediets that the resonanees 
are surrounded by thiek stoehastie layers, while in the density map the resonanee layers 
appear to be stable (and thus lead to higher mean densities). This eomes from a failure of 
frequeney determination: Near resonanees it ean take extremely long integration periods to 
ensure aeeurate reeovery of the frequeneies due to aliasing (Laskar 2003). However, at the 
edges of the resonanee layers—where we expeet to ffnd stoehastie layers and do see weakly 
ehaotie motion—the frequeney diffusion time is a robust indieator of ehaos. 

We eonelude from these experiments that the degree of ehaos is an indieator that en¬ 
sembles of orbits may mix faster than predieted from regular phase-mixing, however the 
Lyapunov time and frequeney diffusion time are not good predietors for the timeseale over 
whieh this mixing will oeeur. To understand this diserepaney, we next explore why this 
oeeurs (Seetion 4.3, in the eontext of Seetion 2.4). 
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4.3. Part III: Short-time frequency evolution 

Why does chaos manifest itself after just 64 orbital periods around orbits with predicted 
chaotic timescales larger than thousands of orbital periods? Both the Lyapunov exponent 
and the frequency diffusion rate measure mean, long-term rates of chaotic diffusion. If 
a weakly chaotic orbit is conffned (by other nearby, non-overlapping resonances or stable 
regions), the mean drift of an orbit in frequency space may be small if computed over 
timescales long compared to the orbital time but short compared to the Arnold diffusion time, 
however small-scale variation of the frequencies may occur over much shorter timescales. 

As a demonstration of the small-scale frequency modulation, we consider again the four 
orbits of Section 4.2 (initial conditions are listed in Table 2). To resolve the short-time 
behavior of the frequency diffusion (corresponding to motion across or around resonance 
layers) we compute the frequencies in a series of rolling windows along each orbit. We use a 
window with a width equal to 64 orbital periods and shift the window by an orbital period 
between each calculation of the most signiffcant frequencies. Figure 11 shows plots of the 
frequencies of each orbit computed in each window of time—plotted in projection are the 
percent deviations of the frequency from the initial value. The difference between the first 
window (blue +) and the last window (red x) is 128 orbital periods. For the two regular 
orbits, (A) and (B), the frequencies vary only from numerical issues when recovering the 
frequencies (^1%) and thus the initial and final value markers overlap. The frequencies of 
the weakly chaotic orbit (C) evolve quickly but are bounded to a small volume with fractional 
size ^1% (presumably by nearby resonant surfaces); this is larger than the frequency spread 
in globular cluster tidal debris (0.1-0.5%). The frequencies of the strongly chaotic orbit (D) 
also evolve quickly but fill a larger volume. 

We see now where the discrepancy between chaotic timescales and the observable effects 
of chaos arise in tidal debris: even if the large-scale diffusion of frequencies is slow, an 
orbit may explore a region of frequency space over much shorter times. The frequency 
diffusion time (Equation 12) is an estimate of the time over which the mean value of the 
frequencies evolves, however the variance of the frequencies over short times is dynamically 
relevant for small ensembles. This small-scale variability is insignificant for the evolution of 
global structure in, for example, an elliptical galaxy or in erasing substructure in the Solar 
neighborhood (e.g., Mafffone et al. 2015), but is signiffant for the morphological evolution of 
tidal debris with small spreads in frequencies. 

We therefore expect a small ensemble of orbits in frequency space to expand over short 
times around even weakly chaotic parent orbits and the debris should therefore appear 
dynamically hotter in real-space. The effect is especially signiffcant if the chaotic evolution of 
the frequencies occurs orthogonal to the largest eigenvector of the distribution of frequencies 


or the local Hessian of the Hamiltonian. We illustrate this phenomenon by computing 
the frequencies for each orbit in small ensembles of orbits around each of the four orbits 
(A,B,C,D). 

We generate orbit ensembles with 128 orbits around the four orbits and integrate for 
two consecutive windows of 128 orbital periods. Figures 12 and 13 show the distribution of 
frequencies for all orbits in the ensembles for the four orbits in each of the two consecutive 
windows. Around the near-resonant orbit (A), the frequency distribution is nearly one¬ 
dimensional, which explains the thin, one-dimensional morphology of the ensemble in the 
left column of Figure 7. The non-resonant orbit ensemble (B) is ‘thicker’—the ratio of 
the two largest eigenvalues of this distribution is larger—and therefore in real space, the 
ensemble begins to spread over two dimensions, as is also evident in the hnal ensemble 
debris morphology in the middle-left panel of Figure 7. Non-resonant orbits in axisymmetric 
or triaxial potentials will generically have frequency distributions that have more comparable 
spreads in each direction and the conhguration-space density of typical tidal debris will 
therefore decrease faster in such potentials (Helmi & White 1999). The weakly chaotic 
orbit ensemble (C) is clearly more diffuse than the two regular ensembles: As is evident in 
Figure 11, the small-scale chaotic evolution of the frequencies—though bounded—increases 
the variance of the frequency distribution. By comparing Figure 12 and 13, it appears that 
the frequency distribution gets more diffuse until it uniformly fills the allowed volume. The 
strongly chaotic orbit ensemble (D) already has a large spread in Figure 12, and in Figure 13 
many of the frequencies have fractional frequency differences relative to the parent orbit 
?:i5-10% (off of the plot). 
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Fig. 10.— The same grid of orbits as shown in Figure 3, but now the greyseale intensity is set 
by the mean density of a of small ensembles of orbits after integration for 64 orbital periods. 
The orbit ensembles are generated around eaeh parent orbit in the grid (Seetion 4.2) and the 
mean density is a kernel density estimation (KDE) with an adaptive Epaneehnikov kernel. 
Given the extremely long ehaotie timeseales of the weakly ehaotie orbits in the upper half 
of the grid, it is surprising that any of the high-order resonant strueture is visible in this 
map. Eour orbits are pointed out with arrows—from top to bottom, these are the strongly 
ehaotie (D), near-resonant (A), weakly ehaotie (C), and non-resonant (B) orbits of Table 2 
(see Seetion 4.2). 
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Fig. 11.— Evolution of the fundamental frequeneies eomputed over a total of 256 orbital 
periods for the four orbits, A, B, C, D. Plotted are the per eent deviations of the frequeneies 
eomputed in window relative to the value in the initial window—that is, if j is the index of 
a given time window, SQij = (Qij — The initial value is shown as a blue plus 

sign and the hnal value is shown as a red x. The frequeneies are eomputed with a window 
width of psl28 orbital periods, and the window is shifted by one orbital period between eaeh 
eomputation (eaeh grey point represents the fundamental frequeneies eomputed in a single 
window). For the regular orbits (A,B), the fraetional variation is around 10“® and thus 
all points overlap on this seale. The weakly ehaotie orbit (C) displays frequeney variations 
eomparable to the spread in frequeneies in globular-eluster-like tidal debris. The frequeney 
variations of the strongly ehaotie orbit (D) have a larger eharaeteristie spread. 























- 31 - 


•cS> 

CN 

G 


2 %- 

0 %- 

- 2 %- 


•!S> 

CO 

G 


2 %- 

0 %- 

- 2 %- 


A 

~r“ 


/ 


-2% 0% 2% 

SQ,, 


0-128 T 
B C 


/ 




- 2 % 0 % 2 % - 2 % 0 % 2 % 


D 

1 -r“ 


. 






<5f*u 


- 2 % 0 % 2 % 

5f!i, 


Fig. 12.— Distributions of the fundamental frequencies for all orbits in 128-orbit ensembles 
generated around the four orbits, A, B, C, D. Plotted are the per cent deviations of the 
frequencies of each ensemble orbit relative to the frequencies of the parent orbit—that is, 
if i is the index of a given orbit and 2 = 0 is the parent orbit, = (Oi^^ — Oi^o)/^i,o- 

All orbits are integrated for 128 orbital periods to compute the frequencies. The near¬ 
resonant ensemble frequency distribution (A) is nearly ID, and thus the ensemble appears 
one-dimensional in conhguration-space (Figure 7). The two largest eigenvalues of the non¬ 
resonant ensemble frequency distribution (B) are closer, and thus the ensemble spreads 
quicker, leading to a two-dimensional spread of debris in conhguration-space. Around both 
the weakly chaotic orbit (C) and strongly chaotic orbit (D), chaotic diffusion increases the 
spread of the debris signiffcantly, leading to faster 3D spreading in conffguration-space. 
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Fig. 13.— Same as Figure 12 but after integrating the same orbits for another 128 orbital 
periods. The ensemble frequeney distributions around the regular orbits (A,B) appear nearly 
identieal, whereas both the weakly and strongly ehaotie ensemble frequeney distributions 
spread signiheantly. 
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5. Discussion: limitations and future work 

We have shown that the Lyapunov and frequeney diffusion times are indieators of ehaos 
and that the frequeney diffusion time resolves the detailed resonant strueture of gravitational 
potentials, but the timeseales predieted do not eapture the importanee of ehaos for the 
density evolution of ensembles of orbits meant to mimie tidal debris. We have shown that 
small-seale but fast ehaotie diffusion of frequeneies ean explain this enhaneed mixing. In the 
sub-seetions below, we diseuss a few important limitations that remain for exploration in 
future work. 
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5.1. The progenitor mass scale 

We have only eonsidered low-mass progenitor systems sueh as globular elusters beeause 
the intrinsie spreads in fundamental frequeneies are small (0.1-0.5%). Small ehanges to 
the frequeneies of the orbits of tidal debris stripped from these progenitors due to ehaotie 
proeesses will therefore eause observable ehanges to the real-spaee morphology of the debris. 
For more massive progenitor systems, the typieal size and veloeity dispersion of the debris 
will be larger and thus the debris morphology will be less sensitive to small ehanges in orbital 
properties. The mass-seale of the debris that will display enhaneed density evolution depends 
on the magnitude of weak ehaos, whieh depends on the orbit strueture of a given potential. 
In potentials with more signiheant ehaos, debris disrupted from more massive progenitor 
systems may also display ‘stream-fanning.’ However, sinee the seale of the debris is larger it 
is more likely that some orbits will exist in ehaotie regions. 


5.2. Potential choice 

The potential eonsidered in this work is ‘unrealistie’ for the total potential of a Milky- 
Way-like galaxy in that it is statie, smooth, and does not eontain baryonie eomponents. 
Simulations of forming galaetie disks in eosmologieal dark-matter haloes have shown that 
baryonie feedbaek and relaxation ean signiheantly ehange the inner distribution of dark mat¬ 
ter and either make the potential more spherieal or oblate (Dubinski 1994; Kazantzidis et al. 
2004; Bryan et al. 2013; Butsky et al. 2015). However, the signiheanee of baryonie relax¬ 
ation or of time-dependenee, triaxiality, and substrueture on shaping the matter distribution 
within the Milky Way is largely unknown. Here we briefly summarize future direetions for 
potential models: 

• Baryonie components: (Debattista et al. 2008, DOS) and (Valluri et al. 2010, VIO) 
studied the orbit evolution indueed by growing a baryonie disk in dark-matter haloes 
with various shapes and orientations (e.g., prolate, triaxial). In general, the authors 
hnd that the growth of a disk slowly deforms the orbits of mass traeers (e.g., dark 
matter partieles) and preferentially populates tube or other ‘round’ orbits (i.e. even 
the ehaotie and box orbits hll a fairly spherieal or oblate volume). Consequently, the 
inner shape of the potential beeomes more oblate or spherieal. If the inner haloes of 
galaxies are indeed elose to spherieal or oblate, the majority of orbits will be regular 
and ehaos will be less important, however this is far from eonelusive. We have found 
from simple experiments in superposing potential eomponents that transition regions 
between potential shapes ean signiheantly enhanee the amount of ehaos. Though su- 
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perposing potentials is not realistic, it at least suggests that a more careful exploration 
of potential configurations is required to understand how complex, radius-dependent 
potential forms affect the amount and signiffcance of chaos in galactic haloes. 

• Time dependence: Galaxies are certainly not static systems. To first order, galaxies 
grow in mass—for example, the spherically averaged mass profile of dark-matter haloes 
evolves fairly predictably in cosmological simulations (Wechsler et al. 2002; Buist & 
Helmi 2014) after some initial period of stochastic mass growth. The Milky Way has 
probably had a fairly calm accretion history over the last 6 Gyr and therefore the 
mass growth may be similar to that seen in simulations. This steady growth most 
likely does not alter the global structure or shape of the potential. However, we also 
know from simulations that figure rotation, baryonic feedback, and the accretion and 
phase-mixing of subhaloes do perturb the global state of simulated haloes. Deibel 
et al. (2011) showed that by adopting pattern speeds comparable to those found in 
cosmological simulations, figure rotation generally acts to destabilize orbits (rather 
than stabilize chaotic orbits in the equivalent static potential) and the resulting orbit 
structure is that most regular orbits are associated with resonances. The presence of 
and response to the Large and Small Magellanic Glouds may also introduce significant 
(time-dependent) perturbations to the global potential of the Milky Way (e.g., Besla 
et al. 2010; Gomez et al. 2015). In future work we will explore the effect of these 
time-dependent processes on the chaotic dispersal of tidal streams using live potentials 
from cosmological N-body simulations. 

• Substructure: Gosmological simulations predict that dark-matter haloes are filled with 
substructure in the form of dark matter subhaloes. If they exist, these subhaloes may 
account for up to ^1-10% of the mass of the dark matter (e.g., Diemand et al. 2007) 
and therefore may contribute significantly to and orbit in the large-scale potential of 
any galaxy. Gravitational scattering due to subhalo interactions has been studied, 
however in the haloes of galaxies where dynamical times are long, the scattering cross- 
section for strong encounters is small. Instead, the collective effect of the subhaloes may 
instead act as a noise term in the Hamiltonian of any halo orbit. This subhalo-induced 
heating—which depends on the mass spectrum and distribution of subhaloes—may also 
act to simply increase the magnitude of chaos along orbits, and destabilize sufficiently 
non-resonant orbits (see, e.g., Kandrup et al. 2000; Siegal-Gaskins & Valluri 2008). 
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5.3. Stream modeling 

Tidal stream modeling is one of the most promising ways to eonstrain the 3D mass 
distribution around the Milky Way at distanees of 10s to 100s of kpe. Methods that use tidal 
streams to infer properties of the Galaetie potential typieally operate by eonstrueting models 
of the debris distribution using either the present-day phase-spaee density (e.g., Varghese 
et al. 2011; Kiipper et ah 2012, 2015; Gibbons et ah 2014), time-of-disruption phase-spaee 
density (Priee-Whelan & Johnston 2013; Priee-Whelan et al. 2014), or the density in angle- 
aetion eoordinates (Sanders 2014; Bovy 2014). All of these methods may fail or produee 
uninterpretable results if modeling globular-eluster streams on mildly ehaotie orbits. For 
eaeh of these methods, it is important to understand the failures and biases introdueed by 
ignoring ehaotie orbit evolution. 


5.4. Ophiuchus stream 

The Ophiuehus stream (Bernard et al. 2014; Sesar et al. 2015) appears to be a thin, 
short tidal stream (deprojeeted length ^1.5-2 kpe) near the Galaetie bulge with no apparent 
progenitor. At Galaetoeentrie R ^ 1 kpe, ^ ^ 5 kpe, the orbits of the stream stars likely 
pass through the MW disk, feel the triaxiality of the Galaetie bar (e.g., Wegg & Gerhard 
2013; Wegg et al. 2015), and have short orbital periods (relative to streams in the halo). It 
is possible that the observed debris is the last remnants of the reeently disrupted progenitor 
system (Sesar et al. 2015), however if a signiheant number of stars were stripped on previous 
perieentrie passages, this older debris may be ‘fanned’ and still near the observed portion 
of the stream. The fanned debris would have signiheantly lower surfaee brightness and thus 
would be more difheult to deteet. The deteetion of this low-surfaee-brightness eomponent 
would open up the possibility that enhaneed density evolution due to ehaos is a dynamieally 
relevant proeess for thin streams in real galaxies, though Garlberg (2015) has shown that it is 
also possible to get short, high density segments of streams from debris formed on eeeentrie 
orbits. 


6. Summary and Conclusions 

We have eonsidered here a simple triaxial gravitational potential ehosen to mimie the 
median properties of dark-matter haloes formed in dark-matter-only simulations (or the 
large-seale properties of haloes formed in simulations with baryonie effeets). We have nu- 
merieally eomputed the magnitude of ehaos for a large grid of iso-energy orbits in this 
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potential using two independent methods that have been used extensively to elassify and 
eharaeterize ehaotie orbits: 1) the Lyapunov exponent and 2) the frequeney diffusion rate. 
From eaeh of these indieators, we eompute a timeseale over whieh ehaos is likely to be impor¬ 
tant and ffnd that the majority of orbits have ehaotie timeseales greater than 100s of orbital 
periods, however with the frequeney diffusion rate we are still able to resolve weak ehaos. 
We then study the density evolution of small ensembles of orbits generated around eaeh 
orbit in the grid used for the previous experiment and ffnd that along some orbits elassiffed 
as weakly ehaotie—with ehaotie timeseales of 100s of orbital periods—the orbit ensembles 
display enhaneed density evolution and reaeh a lower overall density faster than orbit en¬ 
sembles around nearby regular orbits (whieh mix due to phase-mixing alone). We explain 
this diserepaney between the predieted ehaotie timeseale and the observed effeets of ehaos on 
tidal debris by eonsidering the nature of ehaotie diffusion: the elassieal ehaos indieators are 
most sensitive to the slow, Arnold diffusion proeess that ean eause large ehanges to orbital 
properties, but small-seale frequeney evolution oeeurs over mueh shorter times as ehaotie 
orbits stoehastieally diffuse aeross the stoehastie layers that surround many resonanees. The 
fundamental frequeneies of a weakly ehaotie orbit therefore vary over a small region (bounded 
by nearby stable resonanees), and when this amplitude is eomparable to or larger than the 
typieal spread in frequeneies of tidal debris from the progenitor, the phase-spaee density of 
the debris will evolve faster to a state of lower density relative to nearby regular orbits. 

Our main results and eonelusions are summarized as follows: 

1. “stream-fanning” of tidal streams on weakly ehaotie orbits (as seen in simulations by 
Pearson et al. 2015) oeeurs due to small-seale ehaotie evolution of the fundamental 
frequeneies of the debris star orbits; 

2. the Lyapunov time and frequeney diffusion time are powerful indieators of ehaos, but 
do not eapture the importanee of small-seale ehaotie diffusion for the density evolution 
of small ensembles of orbits (tidal debris); 

3. tidal debris beeomes diffuse and thus harder to observe on weakly ehaotie orbits when 
the small-seale ehaotie diffusion of the fundamental frequeneies has a seale eomparable 
to the internal spread in frequeneies of the debris; 

4. the details of the enhaneed mixing along weakly ehaotie orbits depends on the resonant 
strueture of the potential; 

5. the eovarianee of the fundamental frequeneies of an orbit over a given time window 
may be a better predietor of the importanee of small-seale ehaotie frequeney diffusion 
on the resulting morphology of tidal debris. 
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Our results provide a elear explanation of how and why the morphology of tidal streams 
alone ean be used to eonstrain the potential of the host galaxy. The longest thin streams are 
most valuable for this effort beeause they have elearly evolved for a long time, but the debris 
remains eompaet. For shorter thin streams it will be hard to deeouple the unknown evolution 
time from enhaneed density evolution from ehaos. The mere existenee of thin tidal streams 
in the halo of the Milky Way either (1) provides useful information about the potential on 
these seales by, e.g., implying a large degree of regularity, or (2) indieates that the thin, long 
streams (e.g., GD-1) are on regular orbits. These are not mutually exelusive—in faet, if the 
streams are on regular orbits, this would be a powerful way to eheek or rule out possible 
potential models by requiring that the progenitor orbits remain regular. 


The authors wish to thank Robyn Sanderson, Dan D’Orazio, David Merritt, and the 
Stream Team for useful eomments and diseussion. The authors additionally thank the referee, 
Walter Dehnen, for insightful eomments that greatly improved this artiele. 

APW is supported by a National Seienee Foundation Graduate Researeh Fellowship 
under Grant No. DGE 11-44155. KVJ and APW were partially supported by the National 
Seienee Foundation under Grant No. AST-1312196 and NASA grant NNX15AK78G. MV is 
supported in part by the University of Miehigan’s Elizabeth Grosby Eund and the Offfee of the 
VP for Researeh and NASA ATP grant NNX15AK79G. AHWK would like to aeknowledge 
support from NASA through Hubble Eellowship grant HST-HE-51323.01-A awarded by the 
Spaee Teleseope Seienee Institute, whieh is operated by the Assoeiation of Universities for 
Researeh in Astronomy, Ine., for NASA, under eontraet NAS 5-26555. 

This researeh made use of Astropy, a eommunity-developed eore Python paekage for 
Astronomy (Astropy Gollaboration et al. 2013). This work additionally relied on Golumbia 
University’s Hotfoot and Yeti eompute elusters, and we aeknowledge the Golumbia HPG 
support staff for assistanee. This work used the Extreme Seienee and Engineering Diseovery 
Environment (XSEDE), whieh is supported by National Seienee Eoundation grant number 
AGI-1053575 (Towns et al. 2014). 


REFERENCES 

Arnold, V. I. 1978, Mathematieal methods of elassieal meehanies 2 

Astropy Gollaboration, Robitaille, T. P., Tollerud, E. J., et al. 2013, A&A, 558, A33 6 

Baffin, J., Kawata, D., Gibson, B. K., et al. 2005, ApJ, 627, L17 3.1 


Belokurov, V., Zucker, D. B., Evans, N. W., et al. 2006, ApJ, 642, L137 1 
Benettin, G., Galgani, L., & Strelcyn, J.-M. 1976, Phys. Rev. A, 14, 2338 A 
Bernard, E. J., Eerguson, A. M. N., Sehlafiy, E. E., et al. 2014, MNRAS, 443, L84 5.4 
Besla, G., Kallivayalil, N., Hernquist, L., et al. 2010, ApJ, 721, L97 5.2 
Belt, R, Eke, V., Erenk, G. S., et al. 2007, MNRAS, 376, 215 3.1 

Binney, J., & Tremaine, S. 2008, Galaetie Dynamies: Seeond Edition (Prineeton University 
Press) 2, B, B, B.l 

Bonaea, A., Geha, M., & Kallivayalil, N. 2012, ApJ, 760, L6 1 
Bovy, J. 2014, ApJ, 795, 95 1, 2.4, 5.3 

Bryan, S. E., Kay, S. T., Duffy, A. R., et al. 2013, MNRAS, 429, 3316 3.1, 5.2 
Buist, H. J. T., & Helmi, A. 2014, A&A, 563, AllO 5.2 

Butsky, I., Maeeid, A. V., Dutton, A. A., et al. 2015, ArXiv e-prints, arXiv:1503.04814 3.1, 
3.1, 5.2 

Garlberg, R. G. 2015, ApJ, 808, 15 5.4 

Ghirikov, B. V. 1960, Journal of Nuelear Energy, 1, 253 2.2 

—. 1979, Phys. Rep., 52, 263 2.3 

de Zeeuw, T. 1985, MNRAS, 216, 273 2.1 

Debattista, V. P., Moore, B., Quinn, T., et al. 2008, ApJ, 681, 1076 3.1, 3.1, 5.2 
Deibel, A. T., Valluri, M., & Merritt, D. 2011, ApJ, 728, 128 5.2 
Diemand, J., Kuhlen, M., & Madau, P. 2007, ApJ, 657, 262 5.2 
Dubinski, J. 1994, ApJ, 431, 617 3.1, 5.2 

Eardal, M. A., Huang, S., k Weinberg, M. D. 2015, MNRAS, 452, 301 1 

Gibbons, S. L. J., Belokurov, V., k Evans, N. W. 2014, MNRAS, 445, 3788 1, 5.3 

Goldstein, H. 1980, Glassieal Meehanies, Addison-Wesley series in physies (Addison-Wesley 
Publishing Gompany) 2, 2.1 


-39- 


Gomez, F. A., Besla, G., Garpintero, D. D., et al. 2015, ApJ, 802, 128 5.2 

Grillmair, G. J. 2006, ApJ, 645, L37 1 

Grillmair, G. J., & Dionatos, O. 2006, ApJ, 643, L17 1 

Hairer, E., Nprsett, S. P., & Wanner, G. 1993, Springer Series in Gomputational Mathemat- 
ies, Vol. 8, Solving Ordinary Differential Equations. I. Nonstiff Problems, 2nd edn. 
(Berlin: Springer-Verlag), xvi+528, a reprinting with eorreetions appeared in 2000. 
3.2 

Helmi, A., k White, S. D. M. 1999, MNRAS, 307, 495 1, 2.4, 4.3 

Hernquist, L., k Ostriker, J. P. 1992, ApJ, 386, 375 4.2, 8 

Hunter, G. 2002, Spaee Sei. Rev., 102, 83 B 

Ibata, R. A., Gilmore, G., k Irwin, M. J. 1994, Nature, 370, 194 1 

Jing, Y. R, k Suto, Y. 2002, ApJ, 574, 538 3.1, 3.1, 1 

Johnston, K. V. 1998, ApJ, 495, 297 4.2 

Kandrup, H. E., & Mahon, M. E. 1994, A&A, 290, 762 2.4 

Kandrup, H. E., Pogorelov, I. V., k Sideris, I. V. 2000, MNRAS, 311, 719 5.2 

Kandrup, H. E., & Siopis, G. 2003, MNRAS, 345, 727 2.4 

Kazantzidis, S., Kravtsov, A. V., Zentner, A. R., et al. 2004, ApJ, 611, L73 3.1, 3.1, 5.2 

Kiipper, A. H. W., Balbinot, E., Bonaea, A., et al. 2015, ApJ, 803, 80 5.3 

Kiipper, A. H. W., Lane, R. R., k Heggie, D. G. 2012, MNRAS, 420, 2700 1, 5.3 

Kuzmin, G. G. 1973, in Dynamies of Galaxies and Star Glusters, ed. T. B. Omarov, 71-75 
2.1 

Laskar, J. 1988, A&A, 198, 341 3.4 

—. 1993, Gelestial Mechanics and Dynamical Astronomy, 56, 191 3.4, 3.4, 4.1, B 

—. 1996, Gelestial Mechanics and Dynamical Astronomy, 64, 115 3.4 

Laskar, J. 1999, in NATO ASI Series, Vol. 533, Hamiltonian Systems with Three or More 
Degrees of Ereedom (Springer Netherlands), 134-150 2, B 


-40- 


Laskar, J. 2003, ArXiv Mathematics e-prints, math/0305364 4.2 
Laskar, J., & Robutel, P. 1993, Nature, 361, 608 3.4 
Law, D. R., & Majewski, S. R. 2010, ApJ, 714, 229 1 
Lee, J., k Suto, Y. 2003, ApJ, 585, 151 3.1, 3.1, 3.1 

Lichtenberg, A. J., k Lieberman, M. A. 1983, Regular and stochastic motion 2.1, 2.2, A 

Lyapunov, A. M. 1992, International Journal of Control, 55, 531 3.3 

Maffione, N. P., Gomez, F. A., Cincotta, P. M., et al. 2015, MNRAS, 453, 2830 4.3 

Merritt, D., k Valluri, M. 1996, ApJ, 471, 82 1, 2.4 

—. 1999, AJ, 118, 1177 2.2 

Moore, B., Governato, F., Quinn, T., Stadel, J., k Lake, G. 1998, ApJ, 499, L5 3.1 
Navarro, J. F., Frenk, C. S., k White, S. D. M. 1996, ApJ, 462, 563 3.1 
Ngan, W., Bozek, B., Carlberg, R. G., et al. 2015, ApJ, 803, 75 1 
Odenkirchen, M., Grebel, E. K., Rockosi, C. M., et al. 2001, ApJ, 548, L165 1 
Papaphilippou, Y., k Laskar, J. 1996, A&A, 307, 427 3.4, B 
—. 1998, A&A, 329, 451 3.4 

Pearson, S., Kiipper, A. H. W., Johnston, K. V., & Price-Whelan, A. M. 2015, ApJ, 799, 28 
1, 2.4, 1 

Pedregosa, F., Varoquaux, G., Gramfort, A., et al. 2011, Journal of Machine Learning 
Research 7 

Price-Whelan, A. M. 2015, SuperFreq, doi:10.5281/zenodo.18787 3.4 

Price-Whelan, A. M., Hogg, D. W., Johnston, K. V., k Hendel, D. 2014, ApJ, 794, 4 4.2, 
5.3 

Price-Whelan, A. M., & Johnston, K. V. 2013, ApJ, 778, L12 5.3 

Prince, P., k Dormand, J. 1981, Journal of Computational and Applied Mathematics, 7, 67 
3.2 


- 41 - 


Romanowsky, A. J., & Kochanek, C. S. 1998, ApJ, 493, 641 3.1 
Sanders, J. L. 2014, MNRAS, 443, 423 1, 5.3 
Sanders, J. L., k Binney, J. 2013, MNRAS, 433, 1813 2.4 
Sesar, B., Bovy, J., Bernard, E. J., et al. 2015, ApJ, 809, 59 5.4 
Siegal-Gaskins, J. M., & Valluri, M. 2008, ApJ, 681, 40 5.2 

Tabor, M. 1989, Chaos and integrability in nonlinear dynamies: an introduetion (New York 
NY: Wiley) A 

Towns, J., Cockerill, T., Dahan, M., et al. 2014, Computing in Science and Engineering, 16 
62 6 

Valluri, M., Debattista, V. P., Quinn, T., & Moore, B. 2010, MNRAS, 403, 525 5.2 

Valluri, M., Debattista, V. P., Quinn, T. R., Roskar, R., k Wadsley, J. 2012, MNRAS, 419 
1951 3.4, 3.4 

Valluri, M., Debattista, V. P., Stinson, C. S., et al. 2013, ApJ, 767, 93 3.1 

Valluri, M., k Merritt, D. 1998, ApJ, 506, 686 2.2, 3.4, 3.4, B 

van Uitert, E., Hoekstra, H., Schrabback, T., et al. 2012, A&A, 545, A71 3.1 

Varghese, A., Ibata, R., k Lewis, C. E. 2011, MNRAS, 417, 198 5.3 

Vera-Ciro, C. A., Sales, L. V., Helmi, A., et al. 2011, MNRAS, 416, 1377 3.1, 3.1 

Vogelsberger, M., White, S. D. M., Helmi, A., k Springel, V. 2008, MNRAS, 385, 236 2.4 

Wechsler, R. H., Bullock, J. S., Primack, J. R., Kravtsov, A. V., k Dekel, A. 2002, ApJ 
568, 52 5.2 

Wegg, C., k Cerhard, O. 2013, MNRAS, 435, 1874 5.4 

Wegg, C., Gerhard, O., k Portail, M. 2015, MNRAS, 450, 4050 5.4 

Wolf, A., Swift, J. B., Swinney, H. L., k Vastano, J. A. 1985, Physica, 285 3.3 

Zemp, M., Diemand, J., Kuhlen, M., et al. 2009, MNRAS, 394, 641 3.1, 3.1 

Zemp, M., Gnedin, O. Y., Gnedin, N. Y., k Kravtsov, A. V. 2011, ApJS, 197, 30 3.1 


s 


-42- 


Zhu, Q., Marinacci, F., Maji, M., et al. 2015, ArXiv e-prints, arXiv: 1506.05537 3.1 


A. Lyapunov exponents 


If one is only interested in eharaeterizing the degree of ehaos, eomputing the full Lya¬ 
punov speetrum for an orbit is often not neeessary. It is usually suffieient to eompute an 
estimate of the maximum Lyapunov exponent by estimating the hnite-time maximum Lya¬ 
punov exponent (FTMLE). Using the dehnition of w from Equation 1, eonsider an orbit 
that is a small deviation away from the parent orbit, w' = w -\- Sw. If the parent orbit is 
ehaotie, the norm of the inhnitesimal deviation should grow exponentially with time with 
some eharaeteristie rate. A, 

||Mi)ll=e"‘||5«’o|| (Al) 

(see, e.g., Liehtenberg & Lieberman 1983; Tabor 1989). Erom this expression, we see that 




(A2) 


where the maximum Lyapunov exponent is the limit as t (X), 


Amax lim A(t). 

t^oo 


(A3) 


Numerieally eomputing this quantity is not trivial beeause (1) obviously the limit to inhnity 
is not possible and (2) the norm of the deviation veetor ||5i(;(t)|| is expeeted to inerease 
exponentially for ehaotie orbits, leading to nonlinear evolution of the deviation and numerieal 
problems. To eireumvent these issues, it is suffieient to instead start a nearby orbit with 
some small initial deviation with norm 5o, integrate for a suffieiently small amount of time, 
T, then renormalize the deviation baek to the initial norm (Benettin et al. 1976). There is 
no general way to determine r exeept to perform eonvergenee tests. The ETMLE after a 
given number of timesteps, A, is then estimated as 


A 



ll<^«;(^i)ll 


(A4) 


where the L are the times at whieh the renormalization oeeurs and the MLE is estimated 
after a very long time to approximate the limit of Equation A3. 


This preprint was prepared with the A AS E-TLX macros v5.2. 
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For most regular orbits, deviations will grow linearly or as a power-law of time. As we 
have seen in Seetion 2, if the orbit is regular, there exists a loeal transformation to aetion- 
angle variables where the angle variables inerease linearly with time, 6i oc 0^- We ean look 
at small variations around the angle-spaee orbit. 


d{50,) dn, 


dt 


dJ, 


SJ, 


These equations are easily integrated: 


Sdi{t) = Sdi{o) 




SJk]t 


It’s elear then that the norm of the deviation veetor grows linearly with time: 

1 1/2 


||5«,W|| = 




n \ ' 0 


1 1/2 


OC t. 


(A5) 


(A6) 


(A7) 

(A8) 

(A9) 


From Equations A2 and A3 it is evident that any deviation veetor that grows as a power 
law with time, will asymptote to 0 from the limit 

In t 

Amax OC lim k-- = 0. (AlO) 

t — i, 

At long times, the numerieally eomputed MLE for a regular orbit should approaeh 0 as 
Eor ehaotie orbits, the divergenee is exponential, and the limit should eonverge to the rate 
of the exponential: the Lyapunov exponent. In praetiee, the MLE is often estimated as 
the mean of A at after the summation diverges from power-law behavior. Eor weakly ehaotie 
orbits, reliable eomputation of the MLE may take integration of thousands of orbital periods. 


B. Super Freq 

The NAEE algorithm operates on numerieally integrated orbital time-series (i.e. posi¬ 
tions and veloeities at a set of equally spaeed times). Complex eombinations of the phase- 
spaee eoordinates of the orbit are Eourier-transformed using a standard EET and the re¬ 
sultant speetrum is then searehed for the strongest frequeney eomponent; this serves as an 
initial guess for the frequeney of a partieular Eourier eomponent. Solving for the frequeney 
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that maximizes the Fourier integral of the FFT with a partieular window funetion allows 
for a more aeeurate determination of the frequeney. It has been shown that this aeeuraey 
eonverges mueh faster than the typieal expeeted for a standard FFT by using a window 
funetion of the form 

OP 

M^(r = f/fint) = (1 +C0S7rrf (Bl) 

where tint is the integration time (Laskar 1999). Onee the strongest frequeney eomponent 
is found, it is subtraeted from the speetrum and the proeess is repeated iteratively. This 
proeedure generates a table of frequeneies for eaeh eomponent of motion whieh must then 
be searehed for the three fundamental frequeneies. The original implementations of NAFF 
have used p = 1 (e.g., Laskar 1993; Valluri & Merritt 1998), but we have ehosen to use p = 4. 
With a higher order window funetion, the eentral peak is broadened, but the amplitudes of 
the side lobes deerease faster (see Hunter 2002); we have found that this additional damping 
of the side lobes allows for more reliable determination of frequeneies from strongly ehaotie 
orbits. 

SuperFreq reeovers the fundamental frequeneies for an orbit faster (with a fewer number 
of terms) when the eoordinates used are ‘elose’ to the angle variables (PL96; Papaphilippou 
& Laskar 1996). PL96 show that a good ehoiee of eoordinates for tube orbits are the Poineare 
sympleetie polar eoordinates, a set of eanonieal eoordinates similar to eylindrieal eoordinates. 
When eomputing the frequeneies for tube orbits, we hrst align the eireulation about the 
z-axis through rotation, transform to Poineare polar eoordinates, then use SuperFreq to 
measure the fundamental frequeneies. We eould equivalently use the Cartesian time series, 
but the eonvergenee of terms is slower (the amplitudes of sueeessive terms deerease slower 
for Cartesian eoordinates). We have tested that our implementation of SuperFreq returns 
the same fundamental frequeneies in either ease for a set of tube orbits. For box orbits, the 
motion is elose to separable in eaeh Cartesian eomponent and we therefore use Cartesian 
eoordinates for estimating the frequeneies for these orbits. 

Figure B.l shows a validation of our frequeney analysis eode in whieh we reproduee 
the frequeney map at a hxed energy of an axisymmetrie, logarithmie potential (Binney & 
Tremaine 2008, pg. 260, Figure 3.45). Plotted are the (Cartesian) frequeney ratios reeovered 
for a grid of iso-energy, box orbits integrated in the potential 

^(x, y,z) = \ In {x^ + {y/0.9f + {z/0.7)^ + O.l) . (B2) 

Following Binney & Tremaine (2008), we generate a grid of orbits on the equipotential 
surfaee ^{x,y,z) = 0.5 (we use a grid with 100000 orbits eompared to their 10000 orbits). 
Eaeh orbit is integrated for ps 40 orbital periods. In this hgure, stable resonanees appear as 
linear over-densities and unstable resonanees appear as linear under-densities or gaps. The 
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Fig. B.l.— A reproduction of Figure 3.45 from (Binney & Tremaine 2008) as a validation 
of our frequency analysis code: Frequency ratios for 100000 isoenergy orbits integrated in a 
triaxial logarithmic potential. Linear features are resonances—stable resonances appear as 
dark lines, unstable resonances appear as linear gaps. 
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regularity of the points in this map reflects the input grid of initial conditions. Points that 
appear to be erratically scattered are chaotic orbits where the frequencies change with time. 



